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Abstract 



We explore the connection between self-organized criticality and phase 

transitions in models with absorbing states. Sandpile models are found to 

exhibit criticality only when a pair of relevant parameters — dissipation e and 

C^ ■ driving field h — are set to their critical values. The critical values of e and h 

C . are both equal to zero. The first result is due to the absence of saturation (no 

i-^ \ bound on energy) in the sandpile model, while the second result is common 

C ' to other absorbing-state transitions. The original definition of the sandpile 

model places it at the point (e = 0, /i = 0"*"): it is critical by definition. 

We argue power-law avalanche distributions are a general feature of models 

with infinitely many absorbing configurations, when they are subject to slow 

driving at the critical point. Our assertions are supported by simulations 

of the sandpile at e = /i = and fixed energy density C, (no drive, periodic 

boundaries), and of the slowly-driven pair contact process. We formulate a 

field theory for the sandpile model, in which the order parameter is coupled 



t"^^ , to a conserved energy density, which plays the role of an effective creation 

rate. 
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I. INTRODUCTION 

Avalanche behavior is common to many physical phenomena, ranging from magnetic 
systems (the Barkhausen effect) |jl| and flux lines in high-Tc superconductors 0, to fluid flow 
through porous media 0], microfracturing processes IQ, earthquakes |^, and lung inflation 
0]. The common feature of all these systems is slow external driving, causing an intermittent, 
widely distributed response. Avalanches come in very different sizes, often distributed as 
a power law. This fact excites the interest of statistical physicists, since power laws imply 
the absence of a characteristic scale, a feature observed close to a critical point. In order to 
describe a critical point, we need only specify a set of critical exponents, whose values are 
determined by general symmetries and conservation laws and do not depend on microscopic 
details of the system. 

Is there a connection between the observed power-law distribution of avalanche sizes and 
critical phenomena? And if so, can we understand the physics of avalanches by applying 
what we know about critical points and universality? A tentative answer to these questions 
was given by Bak, Tang and Wiesenfield (BTW) 0], who proposed that the power laws 
in avalanche statistics are due to a new kind of critical phenomenon, which they called 
self-organized criticality (SOC). In ordinary phase transitions, criticality is attained only 
by fine-tuning certain control parameters (temperature, pressure, etc.) to special values. 
Only close to this critical point is scale invariance observed. BTW suggested that in some 
dynamical systems the critical point is reached automatically, without any fine tuning, thus 
explaining the wide occurrence of power laws in nature. The idea was then exemplified 
by several dynamical models, such as the sandpile 0, and the forest-fire model [§. The 
SOC hypothesis has stimulated an enormous amount of research. If this — as one of its 
originators has concluded — is "how nature works" 0, then the question of "how SOC 
works" becomes all the more urgent. 

The concept of "spontaneous" criticality, as discussed in the SOC literature, presents, 
however, several ambiguities. Several authors have noted that the external driving rate is a 
parameter that has to be fine-tuned to zero in order to observe criticality [10-[T^. On the 



other hand, it would be amazing, without prior knowledge of the critical coupling, to define 
a system like the Ising model so that it is intrinsically at its critical point. SOC appears less 
miraculous if we suppose that there is 'generic scale invariance,' i.e., that criticality obtains 
over a region of parameter space not just a point. But we will argue that as for the Ising 
model, SOC typically exists at a critical point in the relevant parameter space. The identity 
of the parameters has been obscured by the manner in which the models were defined. 
How can a model be critical by definition, when for most statistical mechanics models, we 
don't even know the exact critical point? One way for a system to discover its own critical 
point is through a suitable extremal dynamics, as in invasion percolation; another is that 
we may know the critical parameters a priori, because they are fixed by a symmetry or a 
conservation law, and build these into the definition of the model. 

Recently, a novel mean-field analysis of SOC models was presented |T^, which pointed 



out the similarities between SOC models and models with absorbing states [0,0. (An 
absorbing state is one allowing no further change or activity.) The mean- field theory provides 
a new insight into the origin of SOC, which, in the sandpile model, is essentially the criticality 
of the population of toppling sites. It turns out that SOC corresponds to the onset of non- 



locality in the dynamics of the system. Non-locality, and hence criticality, is obtained by 
fine-tuning the control parameters, precisely as in continuous phase transitions. In this 
paper we focus on the similarities and differences between SOC and models with absorbing 
states. The latter are relatively well-understood: we can identify the order parameter, and 
know the static and dynamic scaling behavior in the neighborhood of the critical point. The 
mean-field and field-theoretic analysis are well-established, and one has some idea how to 
derive these starting from an exact master equation |TB[ . Applying these ideas, we arrive at 
a new understanding of SOC. 

In Sec. II we review the models of interest: contact processes and sandpiles. The for- 
mulation of a general theory of SOC is problematic because of the nonlocal interactions 
implicitly present in these models. For instance, field-theoretic analysis encounters difficul- 
ties related with the singularity of the continuum limit |]T^. Moreover, it is in general not 



possible to treat simultaneously the two timescales of avalanche propagation and external 
driving. Sec. Ill describes how these problems can be solved by a suitable "regularization" 
of the dynamics, which is local in space and time, and presents sandpile criticality as a kind 
of absorbing-state transition. The regularized dynamics readily lends itself to a continuum 
formulation, presented in Sec. IV. This motivates, in Sec. V, a study of sandpiles at the 
critical point, without boundaries or driving, and of the pair contact process subject to a 
slow drive. We present preliminary simulation results of these systems. We summarize our 
new perspective on SOC in Sec. VI. 

II. CONTACT PROCESSES AND SANDPILES 

A. The contact process 

One of the simplest models showing an absorbing-state transition is the contact process 
(CP) 0] (see Fig. |I|a). To each site of a c?— dimensional lattice we assign a binary variable 
(jj = 0,1. (Occupied sites are said to harbor a 'particle.') Occupied sites {a = l) become 
empty (a = 0) at unit rate, while empty sites become occupied with rate wX, where w is 
the fraction of occupied nearest neighbor (NN) sites. The vacuum state (all sites empty), 
is clearly absorbing and is the only stationary state for A < Ac, while for A > Ac there is 
also an active stationary state. (Ac — 3.298 in one dimension.) The order parameter is the 
density pa of active (occupied) sites, which vanishes at the transition as 

Pa ~ (A - Ac)''. (1) 

The simplest (mean-field) description of the CP treats the ai as uncorrelated: 

^ = -(1 - A)P. - Ap^, (2) 

leading to Ac = 1 and j3 = 1. As in equilibrium, we characterize the critical singularities 



by a set of critical exponents [|T9[, such as /5, and z^^, which describes the divergence of the 
correlation length ^ 

e~(A-Ac)-'^-. (3) 



Besides the "thermal" perturbation A = A — Ac, a second relevant field is an external particle 
source h. {h is the rate of "spontaneous" creation at vacant sites.) For A = 0, p^ ~ h^^^'^. 
Other exponents are defined by considering the decay of perturbations to the stationary 
state [T^. Models with a single absorbing state fall generically in the universality class of 



directed percolation (DP), also known as Reggeon field theory p!9|-pl| . 

A more complicated situation arises when many absorbing configurations exist. The 
simplest model to have been studied in detail so far is Jensen's pair contact process (PGP) 
p^ (see Fig. |l]b). In this model, a nearest-neighbor pair of particles may mutually annihilate, 
with probability p, or else, with probability q = 1 — p, create a new particle at a randomly 
chosen NN, provided it is vacant. There are infinitely many absorbing configurations, since 
all that is required is the absence of any NN particle pairs. In one dimension the static critical 
behavior at Qc = 0.9229 is DP-like p^j23[] , but the spreading or avalanche dynamics has 
variable exponents, depending on the particle density in the environment of the seed p4 



A special, "natural" class of absorbing configurations with particle density (f)nat are those 
spontaneously generated by the critical dynamics. DP spreading exponents are recovered 
only if the initial particle density is set to (pnat — 0.242(1) |2^ 

B. The sandpile model 

Sandpile models are cellular automata (CA) with an integer (or in some cases continu- 
ous), variable Zi ("energy"), defined on a (i— dimensional lattice. At each time step an energy 
grain is added to a randomly chosen site, until the energy of a site reaches a threshold Zc- 
When this happens the site relaxes 

Zi >■ Zi Z(^ y^j 

and energy is transferred to the nearest neighbors 

Zj^Zj + y^. (5) 

The relaxation of a site can induce NN sites to relax in turn, if they exceed the threshold 
because of the energy received, and so on. From the moment a site reaches threshold, 
until all sites have again relaxed {zi < Zc, Vi), the addition of energy is suspended. The 
sequence of events during this interval constitutes an avalanche. For conservative models 
the transferred energy equals the energy lost by the relaxing site (J^Vj = ^c), at least on 
average. Usually, dissipation occurs only at the boundary, from which energy can leave the 
system. 

Since the energy input stops during an avalanche, we have, in effect, an infinite timescale 
separation between the toppling dynamics and the external source. Under these condi- 
tions the system reaches a stationary state characterized by avalanches whose sizes s are 
distributed as a power law [[^,p^-|3T|] 



P{s) ~ s-\ (6) 

The model originally introduced by Bak, Tang and Wiesenfeld (BTW) 0| is a discrete 
automaton in which z^ = 2d and yj = 1 (see Fig. |^). 



An interesting variation of the original sandpile is the Manna model ||32| (see Fig. ^. In 



this automaton the critical threshold is 2;c = 2 independent of the dimensionality d, and if a 
relaxation (toppling) takes place, the energy is distributed to two randomly chosen nearest- 
neighbor sites (see Fig. ^). Variations in which part of the energy is kept by the relaxing site 
can also be considered, as well as models in which energy is transferred along a preferred 
direction p8| . 

Finally, sandpile models in which part of the energy is dissipated have been studied |B3 



In continuous-energy models, some fraction of the energy removed from a relaxing site is lost, 
instead of being transferred to one of the neighbors |33|. In a discrete-energy model, such 



as the Manna or BTW sandpiles, one can introduce a parameter e representing the average 
energy dissipated in an elementary relaxation event. The two dissipation mechanisms lead to 
the same effect, namely a characteristic length is introduced into the system and criticality 
is lost. The avalanche size distribution decays as 

P{s) ~ s-^f{s/s,), (7) 

where the cutoff size scales as Sc ~ e^^^^. We can also observe avalanches in the contact 
process, by starting the system with a single particle |^^. The activity may spread over 



many sites before dying out; avalanches are power-law distributed if A is set to its critical 
value. 

At first glance, the BTW sandpile looks quite unlike the CP. One difference is that the 
avalanche dynamics in the sandpile is nonlocal and deterministic. The sandpile model is 
inherently nonlocal because of the implicit time scale separation. A site can receive energy 
only if the system is quescient, i.e., no active sites are present on the lattice. This implies 
that transition rates depend upon the entire set of lattice variables present in the system, 
giving rise to a strongly nonlocal dynamical rule. Given the configuration prior to the 
avalanche, and the location of the newly added particle, deterministic toppling rules govern 
the evolution to the next stable configuration, and this evolution can affect sites anywhere in 
the system. To have any hope of applying the methods used for the CP, we have to assume 
that the deterministic sandpile dynamics can be realized as a limiting case of models with 
local, stochastic dynamics, belonging to the same universality class as the sandpile. (We 
refer to this as "regularizing" the sandpile rules.) The latter hypothesis would need to be 
verified, but seems plausible if the rules respect the same symmetries and conservation laws 
as those of the original model. (We provide an example in Sec. V.) Similar considerations 
apply to 'extremal dynamics,' which requires the action of an omniscient agent to choose 
the next event. Such a dynamics can presumably emerge as a limiting case of local rules in 
which each unit only has information about a finite number of neighbors. 

As originally defined, the sandpile seems to involve no parameters. There is only the 
toppling rule, which, after some time, miraculously yields a critical state. But in devising a 
regularized dynamics, we are forced to include a nonvanishing driving rate by introducing 
the probability h per unit time that a site will receive a grain of energy |T^. (We may fix the 
relaxation rate for active sites at unity.) Energy is distributed homogeneously and the total 
energy flux is given by J^ = hL'^. The parameter h sets the driving timescale or equivalently 
the typical waiting time between different avalanches as t^ ~ 1/h. As h ^ 0, we recover 
the slow driving limit, i.e., during an avalanche the system does not receive energy. This 
formulation of the dynamics has the advantage of being local in space and time. The state 



of a single site depends only on the state of the site itself and its nearest-neighbor sites at 
the previous time step, through a transition probability that is given by the reaction and 
driving rates. 

III. TOWARDS A LOCAL THEORY OF SOC 

After reformulating the sandpile rules as local and stochastic, we can proceed along the 
path followed for nonequilibrium phase transitions. From the master equation we can derive 
mean-field equations that give a qualitative picture of the phenomenon, exploiting several 
analogies with models with absorbing states. The mean- field analysis of regularized sandpiles 
shows that the order parameter is the density pa of active sites (i.e., whose height z > Zc)i 
and that pa is coupled to the densities of 'critical' {z = Zc — 1) and 'stable' {z < Zc — 1) sites 
131 . In mean- field theory, the dependence of the order parameter pa on the parameters h 



and e can be obtained on the back of an envelope. Since energy is conserved in the stationary 
state, the incoming energy flux Jj„ must be balanced by the dissipated energy J out = ^PaL'^- 
From Jin = J out, we obtain 

Pa = -. (8) 

There is no stationary state for h > e; see Fig. ^. The model is critical just in the double 
limit h,e ^ 0, h/e -^ 0, since the zero-field susceptibility x — dpa/dh diverges, implying a 
long-ranged (critical) response function. Critical behavior emerges in the limit of vanishing 
driving field, corresponding to locality-breaking in the sandpile dynamics. The driving and 
dissipation rates are the control parameters of the model; the stationary order-parameter 
naturally vanishes at the critical point. When h = 0, any configuration with p^ = is 
absorbing. Thus there are an infinite number of absorbing configurations for a sandpile. 



just as for the pair contact process. (In close analogy with Ref. [O], in the field theory of 



the PCP the order parameter — the density of nearest-neighbor pairs — is coupled to a 



non-order-parameter field p3 |. 



In absorbing-state transitions, it is very useful to consider the spread of active sites from 

an isolated seed. Following the scaling framework developed by Grassberger and de la Torre 

T9| , we expect that the probability that a small perturbation imposed on an absorbing 



configuration activates s sites scales as 

P{s,e) = s-^g{s/s,{e)), (9) 

where Sc ~ e~^/°" is the cutoff in the avalanche size. The perturbation decays in the stationary 
subcritical state as 

Pa{t)r^t^nt/t,ie)). (10) 

Here t^ denotes the characteristic time which scales as t^. ~ e^^^ii. In this way we have 
translated the avalanche description into the formalism commonly employed to study models 
with absorbing states P^]. 

It is natural to regard sandpile models as having two parameters, e and h, with the 
original models poised, by definition, at the point (0,0"''). It should be evident that e in 
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the sandpile model is a 'temperature-like' variable, playing the same role as Ac — A = —A 
in the CP. For e > we cannot have sustained avalanches; they decay exponentially in this 
subcritical regime. To have self-sustained avalanches, or an active stationary state with h = 
in the sandpile, we would need e < 0, that is to say, the possibility of creating additional 
energy quanta when a site topples. But this immediately raises a new problem: the energy 
will never be lost (except at the boundaries), so in the thermodynamic limit we shall have 
a runaway 'chain reaction' instead of a stationary state for e < 0! The impossibility of a 
stationary state for e < is analogous to the absence of a well-defined free energy in the 
Gaussian model below Tg. Neither model has the saturation effect needed for stability in the 
'low-temperature' phase. In the Gaussian model the stabilizing u(f)^ term is missing from the 
Hamiltonian, while in the sandpile there is nothing to stop energy accumulating. Indeed, 
to do so would mean to lose energy from the system, destroying the conservation law even 
for e = 0. Criticality would then require that e take some negative but a priori unknown 
value. Thus the possibility of not having to tune the system is predicated on the absence of 
saturation, or, equivalently, on having strict energy conservation when e = 0. 

Another example is furnished by removing saturation (i.e., the restriction to at most 
one particle per site) from the CP, resulting in an exactly-soluble birth-and-death branching 
process. In this model, particles disappear at unit rate, and produce offspring at rate A, at 
neighboring sites, whether they are occupied or not. This corresponds to setting 6 = in 
the mean-field equation 

l-(A-l).-^.^. (11) 

The density grows without limit for A > Ac = 1. It is important to note that Ac = 1 not only 
in mean field theory, but in fact for the actual birth-and-death process. Avalanches follow 
power laws, with the survival probability P{t) ~ t~^, for example. Restoring saturation 
(6 > 0) permits the existence of an active stationary state, but at the cost of shifting Ac to 
some larger but a priori unknown value. (The birth-and-death process is free of higher-order 
terms that would renormalize the critical value of the thermal parameter from that given 
by mean-field theory.) As shown in Fig. ^, neither model has a stationary state for negative 
values of the thermal parameter r; this is the main difference from the phase diagram of the 
CP. The sandpile, however, presents a further subtlety: while there is no upper bound on 
energy, the order parameter is subject to saturation, since pa cannot exceed unity! 

The critical point of the birth-and-death process is at A = 1 because this point cor- 
responds to a balance, on average, between births and deaths. Similarly, the sandpile is 
critical at e = because a toppling site sends a particle to each of g neighbors, and each of 
these neighbors is critical with probability 1/g, so the gain and loss terms for the number of 
active sites balance on average. Thus the sandpile and the birth-and-death process have the 
same phase diagram. This does not mean, of course, that the two models share the same 
avalanche dynamics — that of the birth-and-death process is rather trivial. An important 
aspect of the sandpile is that the condition needed for critical avalanches — that a fraction 
1/g of the nearest neighbors of an active site be critical — is established by the transient 
dynamics of the model. The leftovers from preceding avalanches provide the environment 
in which activity is critical. Memory appears to be the crucial feature of SOC models, and 
is due to the presence of a nontrivial threshold for activity; for sandpiles this means that 
^c > 2. (For completeness, we note that a sandpile with Zc = 1 corresponds to a simple 



random walk, with well-known scaling properties. One may think of it as the analog of 
the birth-and-death process, in the family of models obeying strict conservation of particle 
number for e = 0.) The mean- field analysis [|I^ shows that having (on average) a fraction 
1/ g of the nearest neighbors of an active site critical is necessary for having a stationary 
state, in which energy input is balanced by dissipation. That is, the only stationary state 
for a sandpile at (0, 0"*") is a critical state. 

What happens when we impose an activity threshold on the CP? One realization of such 
a threshold corresponds to the PGP. We study the driven PGP in Sec. V. Here we consider 
the field-theoretic description of the PGP |^. As noted above, the order parameter p is 
coupled to a non-order parameter field n representing the density of isolated particles. The 
equations take the form 



dp 



DpV p — ap — bp — wnp + ■ ■ ■ + rjp (12) 



and 



— = DnV'^p + rp- up^ - wnp -\ h ?]„, (13) 

where the noise terms satisfy 

{r],{x,t)iq,{x\t')) =T,^jp{x,t)5{x - x')5{t - t'). (14) 

The field n{x,t) is frozen in regions where p = 0. {li w = 0, Eq. (0) is the minimal field 
theory for the GP |3^.) Now, because of the simple form of the n equation, we can formally 
eliminate this field to obtain 

^ = n VV -dp- bp'' + Vp- wrp{x, t) f dt'p{x, t')e""^*' '''^^"''^ (15) 

ot Jo 



which exhibits a long- memory effect [Q. The nonlocal term turns out to be irrelevant to 
the stationary properties of the active phase: it is exponentially small if the density of active 
sites is different from zero. The situation can be different for spreading from a seed, in which 
case the active sites form only an infinitesimal fraction of the lattice. 

IV. FIELD THEORY OF SANDPILES 

A field theory of sandpiles should parallel that for the PGP in many respects. As noted 
before, a crucial point of the sandpile dynamics is the coupling of the density field pa with 
the background of critical sites pc- Each region devoid of active sites is frozen until such 
a site is generated. The activity spreads and in general alters the configuration before it 
moves away or disappears. The active sites leave a trace of their dynamical history in 
the frozen configurations of critical and stable sites they produce. If new active sites are 
created in the same region at some later time, they will feel the effect of the active sites 
present earlier in the region. This creates a long-range interaction in time and space among 
active sites. The range of this interaction depends on the characteristic timescale of the 
driving, because the fluctuations induced by h destroy the memory effect. Glose to the 
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infinite timescale separation, the characteristic driving timescale diverges and the range of 
the nonlocal interaction extends to the entire system. This picture is valid also for the PCP, 
since in both systems the response function diverges as pa approaches zero and the nonlocal 
term becomes more and more important as p^ — ;> 0. In sandpile models the density of active 
sites is proportional to the external field h and nonlocality is recovered in the limit h ^ 0. 
We turn now to a detailed continuum description of the BTW model. Let Pi(x, t) be the 
density of sites with height i at x. We note that each site is subject to an input of energy 
due to three sources: (1) The external field, h; (2) toppling of active sites at any of the four 
NN's: (4 — e)pa where pa = Z]j>4 Pi and e is the average energy dissipated; (3) a diffusion- 
like contribution: (1 — |)V^Pa- The diffusive term arises because a gradient in pa leads to 
a particle flux: the excess in the mean number of particles arriving at x from the left, over 
those arriving from the right, is jx{x,t) = —(1 — ^d^Pa- The net inflow of particles at x 
is therefore —V ■ j = (1 — |) V^pa- Applying these observations to the mean-field equations 
derived in Ref . [1^ , we can write down the following set of continuum equations: 

-^ = pi+i + (pi_i - Pi) {(4 - e)[pa + ^VVa] + h} + r]l^^ - r]i + r]i.i , < « < 3 (16) 
and 

-^ = -Pi + Pi+4 + (Pi-l - P^){(4 - e)[pa + ^VVa] + h} - TjJ + TjJ^^ - 7]i + 7]i_i , i > A. 

(17) 

(For z = of course, p_i and ri_i are identically zero.) The terms rji represent noise arising 
due to fluctuations in the number of events of a given kind; rji is the contribution associated 
with the reaction i — > z + 1, and rif with toppling: i —* i — A for i > A. Since the number of 
events is Poisson-distributed (approaching a Gaussian in the continuum limit), the variance 
equals the mean, and the noise variance is proportional to the mean rate of the corresponding 
process. Thus we have 

(r/,(x, t)r7,(x', t')) = T6ij6{^ - x')5(t - t')p,{^, t)[(4 - e)pa(x, t) + h] (18) 

and 

(r^Hx, t)r7j(x', t')) = r%5(x - x')(5(t - t')p.(x, t) . (19) 

The noise terms oc V'^pa have been dropped, as they are expected to be irrelevant. 

This set of equations satisfies probability conservation: J2s Ps is a constant, equal to 
unity by normalization. Let C(x, t) = J2sSPs{'^,t) be the local energy density. From Eqs. 
(^) and ( P^D we have 

^ = {l-^)W'pa + h-epa + Vc, (20) 

where 

(r/c(x,t)r/c(x',t')) = r^x-x')(5(t-t')[/^ + epa(x,t)] . (21) 



(Again we have neglected diffusive noise, and have used the fact that in the absence of a 
source and of dissipation, the total energy does not fluctuate.) For h = e = we have simply 



I - VV. (22) 



SO that E = /(i^xC(x, t) is conserved. 



The generalization of Eqs. (|l^) and (|T7| ) to other dimensions, or to other sandpile models 
(e.g.. Manna's), is straightforward, but the analysis of this complicated set of equations is 
problematic. One might try to cut off the hierarchy by simply declaring Pi = for i greater 
than some ic- The choice of cutoff, however, is not obvious, and one would have to add 
suitable correction terms to ensure that energy is conserved when one active site topples 
onto another. (Altering the sandpile rules to forbid such transfers — in effect, constraining 
all sites to have z < Zc — raises an interesting possibility, but one that we shall not pursue 
further here.) As a step toward simplification of the continuum equations, we sum up Eq. 
(0) for i > 4 to obtain 

^ = -Pa + p: + P3{(4 - e)[p, + 1 Wa] + h}- nl + r/: + % , (23) 

where p* = Y.i>% Piy vl — I]j>4 vI: ^^^ Va — ^i>8 vl- Since the density of sites with heights 
> 8 should be negligible, we might ignore the starred terms. Then the active-site density is 



coupled only to ps, identified in Ref. |T3| as the density of critical sites, pc- In that work, 
sites with heights < 3 are considered in a unified manner, as the density of stable sites, p^. 
Eq. (ITB]) shows, however, that the evolution of pc is coupled specifically to p2, not simply 



to ps = Po + Pi + P2- In the mean-field theory |T3|, the quantity u = P2/ Ps is therefore 
introduced. In a spatially homogeneous stationary state, the value of u can be deduced 
from energy conservation. But in the present context, u = m(x, t) is another dynamical 
variable. Thus our attempt to reduce Eqs. ([I6D and (|17|) to a description in terms of three 
basic categories meets with difficulties. 

Rather than pursuing a systematic derivation of a reduced set of equations from Eqs. 
(p^) and (p!7D, we shall use what we have learnt so far, together with the observation that 
in constructing a field theory, a detailed accounting is unimportant, so long as one respects 
the symmetries and conservation laws of the original model. In the present instance, it is 



essential to ensure conservation of energy when e = h = 0. In fact, Eq. (pOD represents this 
explicitly, and shows how the energy density ( is coupled to pa- We therefore retain Eq. 
(pop as one of our basic equations. 

We obtain the other equation by replacing ps in Eq. ( p3|) with /(C) (1 — Pa)'- only non- 
active sites can contribute to the gain term for pa, and they do so at a rate that depends on 
the local energy density. (/(C) plays a role analogous to that of u in the mean-field theory.) 
For small C, far from criticality, one expects the height distribution to be Poissonian, so 
that /(C) oc C^; for large values of C, / will be a decreasing function of C- The case of most 
immediate interest is a system near the critical stationary state, with C — Cc and / ~ pc, 
where Cc and pc represent the average values of energy and the density of critical sites, 
respectively, at the critical point. From the MF solution we have pc = 1/4 for the BTW 
model in the limit /i — > 0, i.e., the density of critical sites is still rather small, and / is an 
increasing function: /(C) = Pc + ^(C ~ Cc) + ■ ■ ■? with A > 0. Presumably, only the linear 
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term need be retained in the vicinity of the critical point. The resulting field theory for the 
regularized BTW sandpile is given by 

^ = -Pa + [Pc + A(C - Cc)](l - Pa){(4 - e)[p, + IvVa] + M + ^a , (24) 

where 

(r7,(x,t)r7„(x',0) = r(5(x-x')(5(t-t')C(x,t)[(4-e)p„(x,t) + /i] , (25) 



together with Eqs. (|20| ) and (|2T|). As in the PCP, our field theory for the sandpile supports 
an infinite number of absorbing configurations: any C(x, t) consistent with p^ = (when 
/i = 0). In both theories, the non-order-parameter field enters the equation for the order 
parameter in the role of an effective creation rate. The crucial difference between the PCP 
and the sandpile is that in the latter case, this auxiliary field is conserved at the critical 
point. 

In a simple mean-field treatment (spatially homogeneous, no noise), we have 

^ = -[^ - «(C - 1)(1 - \)]Pa - [1 + a(C - 1)](1 - \)pI + [1 + a(C - 1)]^, (26) 

and 

|^-.P. + /. (27) 

where we define 4,A{C — Cc) = (^(C — 1) by introducing ^ = C/Cc- For e and h small, 
and h/e « 1, the mean- field equations have the stable stationary solution pa = h/e, 
( = 1 — h/ae. In the case h = 0~^ — the slowly-driven limit /i, e — >■ 0, with h/e -^ — the 
stationary value of ( approaches the critical height Cc- It is easy to recognize then that e 
plays the role of a control parameter, analogous to A in the CP, with the critical point at 
e = 0. 

A different situation is faced when we impose e = h = from the outset, rather than via 



the slow-driving limit. We have from Eg . (pTl) that ( is strictly conserved in this case. The 



average energy density is thus an external parameter that can be freely fixed in the initial 
condition. In this case, in fact, ( is the only control parameter. In the following section we 
present simulations of just such a situation. 

The full analysis of the field theory will be deferred to a future publication. Here we 
simply observe that for e = h = 0, 

C(x, t) = C(x, t = 0) + r rft'VVa(x, t'). (28) 

Jo 

The evolution of the active-site density contains long-memory terms. If C(x, t = 0) = Co(x) ~ 
Cc, then to leading order 

-^ = ^ VVa + a(Co - l)Pa -pI + ^Pa j^ dtV'pa{^, t') + T^a , (29) 

Unlike the PCP, in which the memory terms decay oc exp[—C J dtpa], here the memory 
decays more slowly, via the diffusive relaxation of pa- It is worth noting that even in active 
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regions, fluctuations in the height field ( cannot relax directly; they only do so by inducing 
similar fiuctuations in pa- Relaxation of the latter then redistributes energy along with 
active sites. 

In summary, all of the models discussed so far — sandpiles, the CP, PCP, and the 
birth-and-death process — have a critical point in a space of two relevant parameters, one 
temperature-like (r), the other field-like (h). Criticality requires h = 0. Models such as 
the sandpile and PCP have a nontrivial threshold for activity and therefore exhibit multiple 
absorbing configurations. When such models are run at (r = 0, /i = 0^), then out of 
a range of possible values for one or more non-order-parameter densities (the critical-site 
density in the sandpile, the density of isolated particles in the PCP), the dynamics selects a 
unique value. From this vantage, the fact that certain models are critical by definition is of 
secondary importance. The essential feature is the behavior of a critical system under slow 
drive. We can study a critical, but non-SOC sandpile by setting e = h = 0; conversely, we 
can observe avalanches on all scales in the PCP if we set p = Pc and h = 0^. 

V. SIMULATION RESULTS 

The preceding discussion motivates several new kinds of simulations of the sandpile and 
the PCP. We report some preliminary results in this section. 



A. Sandpiles at e = h = 

In a regularized theory of sandpiles, we need to introduce a dissipation rate e > h to 
realize a stationary state. In the original model, a stationary state is achieved by imposing 
open boundary conditions. While this may be appropriate for modeling processes in which 
stress may only be released at the boundaries of the system, it is an inconvenience theoret- 
ically: it is easier to study criticality in uniform systems; once bulk behavior is understood, 
the effects of various kinds of boundaries can be analyzed. We therefore study a sandpile 
with periodic boundaries. We performed simulations of the stochastic BTW sandpile at 
{e = 0,h = 0). With e = h = 0, the mean-height ( = N/L'^ is strictly conserved; it is an 
additional parameter at our disposal. 

Initial configurations are generated by distributing at random a fixed number A^ of par- 
ticles amongst L'^ lattice sites. (Since the initial configuration is on average spatially homo- 
geneous, all average properties such as densities and correlation functions are translation- 
invariant.) Once all A^ particles have been placed (but not before), the dynamics begins: 
each active site (i.e., having z > Zc = 2d) topples at unit rate. In practice, we maintain 
a list of the current set of Na active sites, choose one at random as the next to topple, 
and update the list following the redistribution of energy to the 2d neighbors. At each 
toppling event, time is incremented by 1/Na — the mean waiting-time to the next event. 
We compute average properties over a set of Ngamp independent trials, each using a distinct 
initial configuration. {Nsamp = 10^ — 10^ depending on the lattice size and the distance from 
criticality.) 

Sustained activity depends upon two factors. First, there must be at least one active 
site in the initial configuration. This condition is trivially satisfied on large lattices, as 
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the probability of having no active sites becomes exponentially small. (For large L, the 
initial height at a given site is essentially a Poisson random variable, P„ ~ C"c~''/n!, so the 
probability of having no active sites ~ (1 — P2d)^ •) The second requirement is that there 
should be on average at least one critical site amongst the nearest-neighbors of an active 
site. One expects the latter condition to depend sensitively on C,., raising the possibility of a 
phase transition as we vary this parameter. 

In one dimension, not surprisingly, we observe a rather simple behavior. For N < L, all 
trials die out rapidly, so that the only stationary state is the vacuum. For N > L, on the 
other hand, virtually all trials survive indefinitely. (We verified this up to L = 1000. In 
some instances the system evolves to a configuration of the form ...11112011111... in which 
the active site must forever circulate.) Thus we see a first-order transition at C = 1; the 
stationary active-site p^ density jumps from zero to about 0.15. 

In two dimensions the non-driven sandpile exhibits a critical point. Fig. |^ shows that the 
active-site density in surviving trials exhibits a non-monotonic approach to its stationary 
value. By performing studies of this kind, always being careful to check that the system 
has reached a stationary state, we determined p^(C, L) for a range of ( values and for L = 
20, 40, 80 and 160. In Fig. ^ we see that J)^ appears to increase continuously from zero at 
a critical value of (. To fix (c we study the dependence of p^ on L, as it should follow a 
nontrivial power-law (p^ ~ j^-P/'^i- jn the usual notation), only at the critical point. Fig. |^ 
shows 'p^ fall owing a power-law for C, = 2.125, but clearly not for 2.124 or 2.126, allowing us 
to conclude that (c = 2.1250(5) for the two-dimensional sandpile. Indeed, this value for the 
mean height is in perfect agreement with the exact result ( = 2.1248... derived by Priezzhev 
for the driven sandpile |^. We also verify that at C = Cc, each active site has, on average. 



one critical nearest neighbor. The overall density of critical sites is p^ = 0.434, again in 
agreement with driven sandpile simulations [^ (At the critical point, about 10% of critical 
sites have heights in excess of 4.) 

Having located the critical point, we can examine the critical scaling of various quantities. 
Fig. 1^ shows a clear power-law dependence of the active-site density on the distance from 
the critical point: p^ ~ (C ~ Cc)^ with /3 = 0.59(1). The dependence of on JJ^i^^c, L) on 
system size yields /3/z/^ = 0.67(1). (Figures in parentheses denote two standard deviations 
in a least-squares linear fit.) We also monitored P{t), the fraction of surviving trials at 
time t. (Approximately half of the trials appear to survive indefinitely at (c-) Associated 
with the (approximately exponential) approach of P{t) to its limit is a relaxation time, r. 
We find that r has a power-law dependence on L at the critical point: r ~ L~'^^^''^^, with 
^\\/'^± = 1.86(8). For comparison, we note the values for DP in 2-1-1 dimensions: /3 ~ 0.58, 
/3/z/^ ~ 0.80, and 'y\\/iy± — 1.76. The similarity in /3 values is curious, but the differences 
in the other ratios indicate that the sandpile is not in the DP universality class. (This is 
as expected, given the differences between the sandpile and the CP discussed in Sec. III.) 
Studies of correlation functions, that will allow determination of z/|| and z/^ separately, will 
be reported in a future publication. 

In the simulations just described, we have fixed (, one of the variables that the dynamics 
selects in driven sandpiles with dissipation at the open boundaries. We observe criticality just 
at the value (c observed in the driven case, and other variables such as the critical site density 
assume the same value in the two cases. In effect, we are able to study sandpiles in either of 
two "ensembles," one with fixed energy, the other with this variable adjusted by the system 
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dynamics. Open boundaries, which served, in earher sandpile simulations, as an outlet for 
accumulated energy, are now seen not to be essential for criticality. An independent study 
of energy-constrained sandpiles (again with periodic boundaries), confirms that avalanches 
follow the same power laws as in the original BTW model |^. Finally, we note that our 
observation of criticality — at the same (c as in the BTW model — in a stochastic sandpile 
with fully local rules, supports the expectation voiced in Sec. II, that we can study SOC 
using a regularized dynamics. 

B. Driven Pair Contact Process 

The one-dimensional PGP has a continuous absorbing-state transition at Qc] below this 
value of the creation probability, the system falls into one of an exponentially large (with 
L) number of absorbing configurations, each devoid of NN pairs. In contrast with previous 
studies, here we study a driven PGP. Starting from an empty lattice, we add particles at 
randomly chosen vacant sites, until a NN pair is formed. We then suspend the addition of 
particles, and permit the system dynamics, as described in Sec. II, to operate, until the 
system again falls into an absorbing configuration. We simulate a system of size L = 1000 
with periodic boundary conditions and study the avalanche distributions for different values 
of q, with both parallel and sequential updating. 

We collect statistics on the size and duration of the avalanches for various values of q. 
As illustrated in Fig. H, the avalanche-size distribution P{s) is power-law for some range of 
s, but suffers an exponential cutoff at Sc, which grows as q —>■ qc as 

Sc-^iqc-q)-'^". (30) 

(Note that due to parallel updating, the critical creation rate qc — 0.95 rather than 0.9229 
as found in sequentially-updated simulations.) We see that the slope of the power-law 
distribution is consistent with DP (i.e. r = 1.08). Sequentially- updated simulations (not 
shown) yield r = 1.12 and a = 0.45, very close to the expected DP value of 0.44. In 
addition, we observe that at the critical point, the isolated-particle density approaches its 
natural value, (pnat — 0.2 (parallel updating) (see Fig. H). (Similarly, in the sequentially- 
updated case we observe -^ (pnat — 0.242.) A detailed comparison of avalanche scaling 



under parallel and sequential dynamics will be presented elsewhere [§9 . 

In the slowly-driven PGP, the system dynamics 'self-organizes' the isolated-particle den- 
sity (f) to its natural value, the same as in the non-driven system. This is similar to what 
happens in the sandpile, where the driven system selects the same critical mean height, (c 
that we found in simulations without driving. There is, however, one rather striking dif- 
ference between the models. In the PGP, activity can spread at qc for any cj) in the range 
[0,1/2]. In the sandpile, by contrast, activity cannot spread at all if the critical-site density 
is too low. Each toppling destroys an active site, and at least one of the neighbors must 
take its place if activity is to persist. In the PGP, each particle creation generates at least 
one new pair as well, so the activity has a possibility of surviving even in an empty lattice. 
This suggests that one investigate a modified PGP, in which a pair creates at particle at a 
(vacant) second neighbor, rather than at a NN; in this case new pairs will only be formed if 
(f) is sufficiently large. Other potentially interesting models are a saturation-free version of 
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the PCP, and the PCP in two dimensions, where only two distinct universahty class are pre- 



dicted, namely DP and dynamical percolation |^0[. We defer investigation of these models 
to future work. 



VI. SUMMARY AND PERSPECTIVE 

In this paper we have argued that SOC can be understood as an aspect of multiple 
absorbing-state models under a slow drive. We pointed out the similarities in the phase 
diagrams of the two classes of models (for the sandpile and the birth-and-death process, they 
are identical), and in terms of avalanches and of bulk critical behavior, without boundary 
dissipation. We demonstrated that the sandpile exhibits an absorbing-state transition as 
we vary the mean height, and that the PCP, heretofore studied only as an absorbing-state 
transition, exhibits a power-law avalanche distribution under a slow drive. We also suggested 
several new models to investigate, and derived a field theory of sandpiles. 

Beyond these and other avenues for quantitative investigation, we propose a new view- 
point of SOC itself. What 'goes critical' in sandpiles is pa, the density of active sites. 
The evolution of pa is intertwined with other fields, which are frozen when pa = 0. These 
fields describe an energy density that is strictly conserved at the critical point. In order for 
avalanches to be critical, two conditions are needed. First, the parameters h and e must be 
set to their critical values, i.e., to zero. This is accomplished by the definition of the model, 
rather than by tuning parameters, but seems very similar in principle to criticality in CP- 
like models. The second condition is that the environmental density is such as to support 
avalanches on all scales. Particle conservation plays an essential role in this aspect, with 
the threshold for toppling providing for a certain independence between pa and the overall 
particle density. From this vantage, SOC is an absorbing-state transition riding atop a sub- 
strate that preserves a record of the previous activity. SOC typifies the behavior under slow 
drive, at the critical point of a model with an infinite number of absorbing configurations. 

Finally, we offer a comment on the significance of sandpiles as models or paradigms of 
physical processes. The intention of the remark that the sandpile sits, by definition, at 
the critical point in a two-dimensional parameter space, is not to trivialize it, but rather 
provide insight and access to new conceptual and computational tools. One may argue 
whether there is any point introducing e and h as parameters for the sandpile; we merely 
posit that their discussion seems natural if one wishes to draw an analogy between sandpiles 
and other models with critical absorbing-state transitions. The question 'Why is Nature 
filled with systems that tune themselves to a critical point?" may be replaced with: "Why 
do so many systems share the typical features of conservative, saturation-free dynamics, a 
threshold for activity, and widely separated timescales for external driving, on one hand, 
and above-threshold dynamics on the other." The question of how this facet of Nature works 
remains a deep one. 
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FIGURES 

FIG. 1. a) Transition rates in the one-dimensional contact process. Filled circles denote occu- 
pied sites, open circles, vacant sites; gray sites may be either occupied or vacant, b) Transition 
rates in the one-dimensional pair contact process. 

FIG. 2. The BTW sandpile model. When four grains are accumulated in one lattice site 
(zc = 4), the site relaxes distributing the grains to the neighboring sites. 

FIG. 3. The Manna sandpile model. When two grains are accumulated in one lattice site 
{zc = 2), the site relaxes distributing the grains to two randomly chosen neighbors. 

FIG. 4. Phase diagrams of the sandpile and birth-and-death processes, and of the contact 
process and the PGP. The thermal parameter r corresponds to e in the sandpile, 1 — A in the 
birth-and-death process, and to Ac — A in the contact process, 'nss' denotes a region where no 
stationary state is possible. 

FIG. 5. Evolution of the active-site density pa, and of the density pc of critical sites, in 
the two-dimensional stochastic sandpile at e = /i = 0. System size L = 160; mean height 
C = Cc = 2.125. 

FIG. 6. Stationary active-site density p^ in the two-dimensional stochastic sandpile at 
e = /i = 0, as a function of r = ( — Cc- The inset shows J)^ versus mean height ( on linear 
scales. +: L = 40; x: L = 80; o: L = 160. 

FIG. 7. Stationary active-site density /j^ in the two-dimensional stochastic sandpile at 
e = h = 0, as a function of system size L. o: ^ = 2.124; o: ( = 2.125; □: (" = 2.126. 

FIG. 8. Avalanche-size distribution P{s) in the slowly driven, one-dimensional POP, for various 
values of the creation probability q. The system size is L = 1000 and 10^ avalanches are recorded 
for each curve. 

FIG. 9. The density (j) of isolated particles in the slowly driven, one-dimensional PGP at Qc, 
approaches the natural value. 
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